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SUMMARY 


This  analysis  is  aimed  at  the  near-wall  processes  in  an 
injected,  axisymmetr ic,  viscous  flow.  It  is  a  part  of  an  overall 
study  of  solid  propellant  rocket  instability,  in  which  cold  flow 
simulation  is  evaluated  as  a  tool  to  elucidate  possible 
instability-dr ivinq  mechanisms.  One  such  prominent  mechanism 
seems  to  be  visco-acoustic  coupling,  as  indicated  by  earlier 
detailed  order  of  maqnitude  analysis.  The  major  component  of  the 
overall  study  involves  numerical  simulation  of  the  full  set  of 
coreflow  equations  of  motion  (nonsteady,  axisymmetr ic)  by  a 
modified  MacCormack  integration  technique.  To  clarify  some  of 
the  physical  interactions  inherent  in  the  various  regimes  of  the 
flowfield,  two  (separate)  singular  perturbation  analyses  have 
been  carried  out.  The  head-end  boundary  regime,  and  the  injected 
sidewall  layer,  both  involve  appreciable  viscous  dissipation,  and 
hence  are  characterized  by  predominantly  parabolic  differential 
systems.  The  inverse  square  root  of  the  injection  Reynolds 
number  serves  as  a  small-perturbation  quantity.  The  sidewall 
layer  analysis  yields  a  first  order  axial  pressure  distribution 
(due  to  viscous  effects)  which  correlates  the  available  steady 
state  data  (CSD  experiments,  1982)  very  well.  The  radial 
dependence  of  the  inner  variables  up  to  first  order  is  solved,  so 
that  the  (x,t)  dependence  can  be  described  by  much  simpler 
partial  differential  systems;  inner/outer  matching  was  not 
attempted.  The  finite-difference  coreflow  algorithm  is 
described,  and  a  source  listing  from  one  preliminary  version 
(R0SC0-2)  is  provided.  Convergence  to  steady  state  has  been 
achieved  with  this  program  recently  (by  marching  forward  in  time 
from  imposed  initial  data).  The  modular  structure  of  the 
computer  program  facilitates  addition  of  special  segments  (such 
as  two-equation  turbulence  modeling) . 
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1.  BACKGROUND 

This  is  pact  of  a  study  aimed  at  elucidation  of  the  physical 
mechanisms  capable  of  driving  acoustic  instability  in  solid 
propellant  motors,  particuarly  of  the  type  termed  velocity- 
coupled  instability.  Previous  studies  on  the  coupling  between 
velocity  oscillations  and  the  combustion  process  in  solid 
propellant  motors  have  demonstrated  the  complexity  of  the  overall 
phenomenon,  but  have  not  yet  defined  the  basic  mechanisms  nor  how 
they  operate  under  flow  conditions  prevailing  in  rocket  chambers. 
Critical  literature  review  and  order  of  magnitude  analyses  of 
velocity  coupling  mechanisms  have  been  carried  out,  including 
visco-acoustic  coupling  and  turbulence  combustion  coupling.  The 
major  goal  of  the  study  is  the  analytical  simulation  of  the 
interior  flow  field  within  a  solid  propellant  grain.  The  focus 
is  on  the  Stokes  layer,  with  the  objective  of  investigating  the 
particular  instability  mechanism  of  visco-acoustic  coupling. 
Preliminary  analysis  has  indicated  that  this  mechanism  is  both 
plausible  and  sufficiently  powerful  to  drive  nonlinear 
vibrations?  it  has  been  shown  that  the  frequency-dependent 
surface  heat  feedback  component,  due  to  viscous/acoustic 
coupling,  has  both  phase  amplitude  ranges  which  would  enable 
driving  of  acoustic  vibrations;  its  amplitude  tends  to  increase 
as  the  mean  coreflow  Mach  number  and  the  frequency  become  higher. 
A  comprehensive  analytical  model  of  the  flow  field  within  the 
viscous  wall  layer  region  has  been  derived,  for  an  axisymmetric, 
nonsteadv  flow  field  configuration.  For  a  simulation  of  the  cold 
flow  test  results  gnerated  at  UTC/CSD,  four  conservation 
equations  are  incorporated  for  continuity,  momentum  and  energy. 

The  major  aspect  of  the  near-wall  behavior  from  the  visco¬ 
acoustic  point  of  view,  is  the  laminar  dissipative  processes 
tvpical  to  that  region.  This  analysis  is  focused  on  the  near- 
wall  processes.  Although  the  solultions  derived  are  nonsteady  in 
general,  the  radial  wall-layer  distributions  obtained  could  best 
be  demonstrated  at  steady  state.  For  this  reason,  the  review 
herein  is  limited  to  steady  behavior. 

Culick  [1]  derived  a  solution  to  the  Stokes  stream  function 
equation  or  flow  in  a  pipe  with  injected  sidewalls.  The  flow  is 
rotational,  and  despite  being  inviscid,  could  obtain  a  solution 
for  the  axial  velocity  component  which  satisfied  the  no-slip 
boundary  condition  at  the  wall.  The  solution  which  satisfies  the 
boundary  data,  namely,  u(x*0)*O,  u(r»l)=0,  and  v(r=l)*l,  yields: 


v  =—  Siv>(§-rz)/r^  a  =  irxc«*(f  r1) 


Of  the  general  family  of  solutions  obtainable,  onlv  that  which 
allows  full  determination  of  the  vorticity  (the  azymuthal 
component  alone  remains)  by  the  available  boundary  data,  is 
physically  meaningful;  the  rest  were  therefore  rejected.  The 
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axial  pressure  distribution  obtained  from  the  momentum  equation 
is  parabolic, 

(Po-PVfv1  =  fx3- 

This  type  of  injected  flow  field  has  been  investigated  previously 
both  experimentally  and  theoretically.  In  particular  the  early 
theoretical  work  of  Berman  [2],  who  arrived  at  a  power-series 
solution  to  the  perturbation  problem  of  suction  in  a  prismatic, 
porous-walled  channel,  with  the  suction  Reynolds  number  serving 
as  small-perturbation  quantity.  The  analytical  results  of  G. 
Taylor  [3]  and  Wageman  and  Guevara  [4]  more  closely  resemble  the 
cosine  terms  of  Culick  [1];  both  [3,4]  have  carried  out 
experiments  as  well,  and  both  demonstrated  very  good  agreement 
between  the  measured  axial  velocity  profiles  and  the  calculated 
ones.  It  appears  that  Culick  [1]  has  arrived  at  his  results 
independently,  since  no  reference  was  made  to  any  of  the  previous 
works.  In  the  experiments  by  Dunlap,  Willoughby  and  Hermsen  [5], 
the  formulation  derived  by  Culick  [i]  was  used  to  correlate  the 
measured  data,  again  with  considerable  success,  regarding  the 
coreflow  axial  velocity  profile,  that  is,  away  from  the  close 
neighborhood  of  the  wall. 

Other  experiments  by  Olson  and  Eckert  [6]  and  later  by 
Huesman  and  Eckert  [7]  tend  likewise  to  verify  the  validity  of 
this  formulation,  in  particular  regarding  the  radial  velocity 
profile,  which  indeed  exhibits  a  peak  near  the  porous  surface 
[6],  as  well  as  the  axial  pressure  distribution  (the  latter  shown 
as  a  linear  correltaion  between  the  friction  coefficient,  Cc,  and 
the  inverse  mean  axial  velocity,  which  are  both  proportional  to 
1/x. 


The  recent  (and  ongoing)  experimental  study  by  Brown,  et  al 
[8]  provides  valuable  information  regarding  the  steady  state 
axial  pressure  profile  and  the  axial  velocity  distribution,  as 
well  as  nonsteadv  wall  heat  transfer  (obtained  by  exciting  the 
standing  acoustic  modes  in  the  tube).  Departure  of  the  steady 
state  data  from  the  predictions  of  the  aforementioned  formulation 
by  Culick  [1]  was  attributed  to  possible  transition  to 
turbulence.  As  will  be  shown  in  this  study,  the  pressure  data 
obtained  can  be  simulated  very  well  with  a  first-order  pressure 
perturbation,  arising  from  the  laminar  viscous  wall-layer 
analysis. 

Earlier,  Yagodkin  [9]  reported  an  experimental  cold  flow 
setup,  with  an  injected  porous  pipe.  The  maximal  injection 
Reynolds  number  was  250,  which  is  2-3  orders  of  magnitude  less 
than  that  corresponding  to  actual  internal  rocket  flows.  Hot¬ 
wire  anemometry  was  used  to  obtain  axial  velocity  and  axial 
velocity  fluctuation  vs  axial  and  radial  distance.  Turbulence 
intensitv  seems  to  peak  near  the  surface,  and  decrease  toward  the 
centerline  and  toward  the  pipe  wall.  These  observations  are 
qualitatively  similar  to  those  obtained  later  by  Yamada,  et  al 
[10].  Although  a  transition  region,  at  Reo=100-150,  was 
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soeculated  [9]  to  involve  "large  eddy  structures",  no  such 
evidence  appears  in  the  experimental  data  reported  [9}. 

Further  studies  by  Yagodkin,  with  Varapaev  [11]  and 
Sviridenkov  [12]  are  theoretical,  and  address  the  problem  of 
laminar  stability  of  injected  channel  flows,  i.e.,  transition  to 
turbulence.  Thus,  modified  versions  of  the  Orr-Sommer field 
problem  were  investigated  analytically  [11]  and  numerically  [12]. 
Two  related  laminar  flow  stability  analyses  are  by  Goldshtik,  et 
al  [13]  and  Alekseev,  et  al  [14].  None  of  these  theoretical 
analyses  indicates  the  presence  of  large  turbulent  eddy 
structures  prior  to  a  full  transition  point,  neither  do  they 
obtain  an  origin  of  such  turbulence  on  the  centerline  upstream. 

Recently  Flandro  [15]  has  carried  out  a  theoretical  analysis 
for  a  burning  propellant  in  a  cylindrical  grain,  under  the  effect 
of  incident  acoustic  waves.  A  detailed  formulation  was  derived 
with  a  double  expansion,  in  terms  of  both  inverse  Reynolds  number 
as  well  as  Mach  number  (independent  small  parameters).  A 
nonsteady  premixed  combustion  zone  was  considered  near  the 
propellant  surface;  the  assumption  is  made,  however,  that  flow 
within  the  combustion  zone  is  pure  radial,  i.e.,  zero  axial 
component  to  all  orders.  Thus  it  could  be  anticipated  that  the 
results  resemble  (regarding  nonsteady  combustion  behavior)  those 
of  T'ien  [16],  and  there  seems  to  be  only  small  differences 
between  the  response  to  tangent  and  to  perpendicular  wave 
incedence.  The  problem  is  finally  solved  numerically,  and 
details  of  the  inner/outer  matching  process  were  not  given. 
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2.  ANALYTICAL  MODEL  OF  THE  COREPLOW 
2 . 0  INTRODUCTION 

In  this  chapter,  the  equations  of  motion  pertaining  to  the 
core-flow  simulation  are  presented,  for  an  axisymmetric  flow 
field.  Turbulence  and  combustion  are  precluded  from  the  present 
formulation,  for  reasons  discussed  earlier.  Other  than  these 
simplifications,  the  full  compressible,  noneteadv,  viscous 
equations  of  motion  are  considered,  with  all  the  dissipative 
terms  included.  A  schematic  of  the  simulated  flowfield  with  the 
various  regions  of  interest  is  shown  in  Fiqure  1. 

Treatment  is  divided  into  three  subsections,  in  which  the 
coreflow  region,  the  head-end  closure  region,  and  the  porous, 
injection  sidewall  region,  are  discussed  in  detail.  The  latter 
two  parts  represent  singular  perturbation  analyses  of  the 
boundary-layer  type,  which  are  incorporated  for  generation  of 
boundary  data  within  the  coreflow  solution  procedure.  Important 
physical  insiqhts  are  obtained  regarding  the  behavior  of  the 
system  at  relatively  large  injection  Reynolds  numbers,  coupled 
with  low  injected  Mach  numbers,  such  that 

Mq2  ~  O(l/Reo) 

The  numerical  algorithm  developed  for  solution  of  the 
coreflow  differential  system  is  a  modified  MacCormack  scheme. 

Its  finite  differencing  details  §re  discussed  in  the  next 
chapter,  along  with  a  current  listing  of  the  Fortran  code. 
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2.1  THE  COREFLOW  FORMULATION 


The  objective  is  to  simulate  the  cold-flow  experiments  of 
Dr.  Brown  at  UTC/CSD,  which  utilize  cylindrical  geometry.  For 
this  purpose,  an  axisymmetric  formulation  was  derived,  to 
describe  a  nonsteady,  compressible,  viscous  flowfield.  For  the 
coreflow  region,  with  typical  injection  Reynolds  numbers  of  order 
1000  and  larger,  we  assumed  constant  and  uniform  thermophysical 
properties.  As  mentioned  earlier,  combustion  (or  chemical 
change)  and  turbulence  are  precluded  at  the  present  stage. 

The  five  equations  of  motion,  for  continuity,  radial 
momentum,  axial  momentum  and  energy  are  presented  in  differential 
form.  A  caloric  eouation  of  state  (pertaining  to  perfect  gas) 
completes  the  model  to  form  closure  of  the  dependent  variables. 

The  following  dimensionless  independent  variables  are 
introduced,  based  on  the  two  physical  scales  of  reference, 
inner  chamber  radius,  R  *,  and  reference  injection  velocity, 

vo*: 


r  =  r 


‘/v 


,  X  =  X*/R  *  | 


t  =  t*/tQ* 


(2.1) 


where  tQ*  =  R0*/v0* 

The  dependent  variables  are: 

v  =  V*/V0*,  u  =  u*/vQ*,  hQ  =  h*/h0* 


(2.2) 


and 


P  =  P*/Pn* 


(2.3) 


In  the  last  equations,  the  properties  used  for  non- 
dimensionalization  are  the  reference  (injected)  density,  $>0*,  and 
the  reference  chamber  pressure,  pQ*;  the  corresponding  thermal 
enthalpy,  hQ*,  is  calculated  from  the  caloric  equation  of  state. 


TT-{ 


o 


(2.4) 


where  =Cp/Cv  is  the  specific  heat  ratio,  is  considered  as 
constant.  The  reference  speed  of  sound  is 


a:  =  ( -*?:/?:  =  v?<>ow 


.* 

o 


(2.5) 


The  corresponding  injection  Mach  number  is 

M  =  v  *  /a  * 

w n  t 


(2.6) 


The  reference  (injection)  Reynolds  number  and  Prandtl  number  are, 
respectively. 


Reo  =p0*v.*e*//u* 

%  =  I iScf/*.* 

Recall  that  the  viscosity,  thermal  conductivity  and  isobar ic 
specific  heat  are  all  uniform  and  constant  within  the  present 
cold-flow  simulation.  These  idealizations  are  incorporated 
merely  for  convenience,  in  allowing  clear  identification  of 
physical  interactions  within  the  coreflow,  at  low  (axial)  Mach 
numbers?  sharp  pressure  and  temperature  variation  are  obviously 
precluded.  These  simplifying  assumptions  are  in  no  way  essential 
to  the  solution;  variable  viscosity,  thermal  conductivity,  etc., 
can  be  radily  incorporated  in  the  numerical  solution  algorithm. 

The  dimensionless  equations  of  motion  are  as  follows,  for  the 
region 


(2.7) 

(2.8) 


0  <  x  <  1,  0  <  r  <  1,  t>0: 


CONTINUITY: 


‘  L2!-(rfv)-f-  ^  =S=o 


2 v  +  r  -or’ 


(2.9) 


RADIAL  MOMENTUM: 


•at 


(2.10) 


AXIAL  MOMENTUM: 


7ft  r 


Sa ) =  £ 


(2.11) 


THERMAL  ENTHALPY: 


L  — ■  (t  p  phv)+  h  it)  = 


(2.12) 


2ft  r  dr 
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where  the  right-hand  side  (source)  terms  are  defined: 


4/3  J_/'2y  _y.'\  .  < 

trt*r  rJ+ee0{  ^+*>  vm. 


\ 

3  •^r2-/ 


(A:  f  m 

Trr^T  'arj 


■t- 


-L  CdJ£+f2Ju  + 

&ec,  fT1-  3  a** 


'dTdxJ 


(2.13) 


(2.14) 


(2.15) 

The  parameters,*^,  Pr,  MQ2,  Reo  are  all  constants. 

The  differential  system  for  the  coreflow  can  be  written 
in  short  notation. 


'VX 


(2.16) 


where  Uk  is  the  non-primitive  dependent  variable  vector,  with  the 
components  defined: 

;>  fl h  (2.17) 

while  the  flux  terms  F^(U),  Gk(U)  depend  only  upon  the  vector  0 
(when  ,  MQ  are  considered  as  constant  parameters: 


(2.18) 


<s>,  =fU-,  G2  =  f  VU  }  64=  f  w.z+ 

^4=r^f>hC>L  (2.19) 

Note  that  incorporation  of  the  pressure  gradient  within  the  S2 
source  term  in  the  radial  momentum  equation,  while  the  axial 
pressure  gradient  is  included  within  the  axial  flux  comonent,  G3, 
is  merely  for  convenience  in  the  solution  process.  In  the 
meantime,  the  viscous  and  thermal  dissipative  terms,  ~0(1/Reo), 
are  expected  to  be  very  small  over  most  of  the  coreflow  domain, 
excluding  the  neighborhood  of  the  walls. 

The  Boundary  Conditions:  The  following  physical  boundary  data 
are  available,  tor  the  cold-flow  simulation: 

(a)  On  the  centerline,  (t,  r=0,  x) : 


< 

II 

0 

V* 

'&A/3T=  iTff  =  0 

(2.20) 

(b)  At  the 

porous  (injected)  surface,  (t,  r*l,  x) : 

v  =  -vQ(x,t),  u=0 ,  h  =  hQ(x,t) 

(2.21) 

(c)  At  the  (nonpermeable  solid)  head-end  closure,  (t,  r,  x=0) : 


v  =  0,  u  «  0,  h  «  hH(r,t)  (2.22) 

The  functions  vQ(x,t),  hg(x,t)  and  hH(r,t)  are  arbitrary  imposed 
(generally  variable)  distributions. 

(d)  The  exit  plane,  defined  by  (t,  r,  x=L) ,  forms  an  entrance 
into  a  short,  convergent  nozzle  section.  This  nozzle  section  is 
treated  separately  from  the  rest  of  the  flow  field;  several 
assumptions  are  incorporated,  as  follows: 

(d.l)  The  throat,  At(t),  is  variable  but  remains 
sonic  at  all  times. 

(d.2)  The  convergent  section  is  short,  and  introduces  no 

dynamic  effect;  it  responds  instantly  to  any  changes, 
and  is  considered  (in  this  sense)  quasi  steady. 

(d.3)  The  sonic  surface  at  the  throat  is  reasonably 

approximated  by  a  plane  encompassing  the  entire 
(circular)  throat  area.  In  other  words,  within  each 
computational  cell,  the  flow  can  be  considered  quasi 
one-dimensional. 

The  last  assumption  is  used  to  facilitate  calculation  of 
the  implicit  functional  relationship 


f (  (u/a) 2 ,  At/A,  .. .]  *  0 


(2.23) 


within  each  computational  cell  in  the  discretized  flowfield  within 
the  nozzle. 

The  foregoing  discussion  has  summarized  the  corelflow 
analytical  model,  including  the  equations  of  motion  and  the 
relevant  boundary  data. 

Simulation  of  the  nonsteady  flow  field,  which  arises  due 
to  perturbation  of  the  exit  nozzle  can  be  performed,  with  the 
initial  data  corresponding  to  steady  state.  As  mentioned 
earlier,  solutions  are  generated  numerically,  by  a  finite- 
difference  algorithm.  Prior  to  the  description  of  the  numerical 
algorithm,  two  special  regions  in  the  flow  are  discussed  in 
detail:  the  head-end  and  the  sidewall  layer,  which  appear  in  the 
following  two  sections. 


2.2  THE  HEAD-END  CLOSURE  LATER 

A  solid,  planar  head-end  closure  is  considered  at  x=0,  as 
shown  in  Fig.  2.  The  flow  in  this  region  is  radially  injected 
inward  (from  the  porous  cylindrical  wall),  then,  near  the 
centerline,  tends  to  turn  toward  the  axial  direction. 

A  boundary  layer  is  formed  near  x=0,  to  connect  the  regular 
flow  regime  (with  apperciable  radial  and  axial  motions),  with  the 
end-wall  where  no-slip  conditions  prevail,  viz.,  v=u=0.  Within 
this  layer  viscous  forces  are  of  importance.  The  dimensionless 
parameter 

0  <  £*  1/  <<  1  (2.24) 

can  serve  as  a  proper  small  perturbation  quantity.  The  layer 
axial  coordinate  is  therefore  stretched, 

yi  =  x/£  (2.25) 

Thus,  the  independent  variable  system  is  transformed  from  (x,  r, 
t)  to  (y r,  t)  in  the  layer.  Further,  the  dependent  variables 
are  now  perturbed,  as 

r-ft+eft  ,  v=Vo«-ev,  s  u.=  u„ (2  26) 

I't  st  Vio  ■*-£!*( 

where  the  dependent  variables  are 

v0(Yi»r»t),  vjJy^rjt),  Uj^yj^r,^  ...  etc.,  all  assumed  to  be  of 
order  unity.  For  convenience,  the  following  abbreviated 
definitions  are  introduced: 

F°  =  j  Vl+5V,,3QlSSJU1+?(Uo  ,2.27, 

The  transformed  axial  derivatives  are  now, 

^/dy.  =  ^  (2.28) 

while  the  (r,t)  variations  remain  equal  to  their  counterparts  in 
the  original  equations  of  motion. 

In  the  remainder  of  this  section  we  will  derive  the 
perturbed  system  of  equations  of  motion  for  the  layer,  and 
collect  hierarchies  of  equal  power  of  &  . 

Obviously,  the  injection  Mach  number  appears  as  an 
additional  parameter  in  the  formulation  (equations  of  momentum 
and  energy).  In  the  flow  fields  of  interest  for  simulation 
herein,  Mq  is  also  very  small;  in  consideration  of  typical 
experiments  at  CSD/UTC  with  air  injection,  we  find 


Mq  -  O(l/Reo)~  EZ 


(2.29) 


which  adequately  represents  a  range  of  cold-flow  conditions. 

This  offers  great  simplification  in  the  analysis,  although  at  the 
cost  of  narrower  range  of  general  application  (considering  the 
relative  freedom  of  the  two  major  flow  parameters,  Reo  and  MQ). 


Therefore  a  oarameter  is  introduced. 


Km  = 


_  '/l2e«  _  S' 


aJ 


:z 

TMo1 

according  to  the  foregoing  considerations. 


O(i) 


(2.30) 


*The  question  of  timescale  is  not  trivial,  since  it  depends 
upon  the  range  of  frequencies  of  interest.  The  following 
reasoning  will  demonstrate  that  for  the  range  of  conditions 
considered  for  the  present  simulation,  the  timescale  can  remain 
the  same  one  as  in  the  coreflow.  The  viscous  layer  thickness  is 


8  v  ~  *o/ V**e o 


where 


®eo  =  v 


(2.31) 

(2.32) 


The  Stokes  Layer  thickness  (for  acoustic  perturbations  with  a 
frequency  f*) 

Sc™  ~  ^/f": 


’STO  ’  */fo 
The  ratio  of  these  two  thickness  scales  is, 


Sy/Ss-ro  ^  \/fo  Ro/V’ 


KO 


’/z 


(2.33) 


(2.34) 


where  SR0  is  the  relevant  (injection)  Strouhal  number.  The  range 
of  Strouhal  numbers  considered  is  SR0  “  0(1),  so  one  need  not 
introduce  any  additional  timescales. 

The  continuity  equation  becomes,  after  using  the  perturbed 
variables,  Eqs.  (2.26)-(2.28) : 

JUsV*  £P.)  +  fc  Vo  ^  e  ^  ^  ^  £<s,+ 


I  .A 


(2.35) 


I. 


From  which  the  following  hierarchy  is  collected! 
ORDER  1/  6  : 


—O 


Gq  =  Gc  (r,t) 


(2.36) 

(2.37) 


ORDER  £.°  (ZEROTH): 


TTc  GT 


lT&_ i 
W)l 


r 


(2.38) 


ORDER  61  (FIRST)  : 


^  3ga,  _  f, 

^  -3y, - r 


Similarly,  the  radial  momentum  equation  yields 


(2.39) 

I;  jb  ft  Mo  £(F,Uc  *-F0  Ui)J  =  —  y=  6 


_  R.U.  ^  V.t+Vit 


(2.40) 


The  following  hierarchy  evolves: 
ORDER  1/  62: 


'dftfyC  =0 


(2.41) 


ORDER  1/  &  s 


XF-u^faj,  +-  =  o 


ORDER  £° (ZEROTH) s 


3k&«0  +  %,(«.*•*  a,e)  =  -^ 


ORDER  61  (FIRST)  : 


+  12 V, 

The  axial  momentum  yields,  after  similar  treatment: 
ORDER  1/  €3:  'dfohiA.  —Q 


ORDER  1/  6  2J 


ORDER  1/fe  : 


'dPobiq(  —  o 

=o 

=  O 


Thus,  in  view  of  Eq.  (2.37) 


uQ  -  uQ  (r , t) 


ORDER  fc°  ZEROTH: 

.  36bVo  *3 


ORDER  eA  (FIRST)  : 


2.5 
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At  this  point,  the  lower-order  anlaysis  results  can  be 
summarized: 

■a^/3r=o  -^34^=0, 


where: 
3  C^oVo 

~W, 


4^  =  f0  h0 . 
£  ^  -n 

Tnr  O' 


(2.51) 


4.1S  -  «.■?§  -0 , 


Based  on  the  foregoing,  along  with  the  boundary  data: 

uo  (Yi=0,r,t)  =0  C 

the  natural  choice  of  solution  for  uQ  is  the  trivial  one, 

uo  (yi»r»fc)  ~  0  (- 

Thus,  GQ(yj_,  r,  t)  =  0.  Further,  from  Eg.  (2.42), 

^P]/c)r  =  0  C 

which  therefore  leaves  only  time-dependent  pressure  within  the 
layer,  up  to  0  (  g,2)  : 

i%  =  fs(.t)  ,  p,  =  r,w  c 

which  implies  also,  for  4>=?k  : 

With  the  foregoing  results  incorporated,  the  energy  equation 
becomes  much  simpler  to  handle;  the  following  hierarchy  is 
obtained: 

ORDER  &°  (ZEROTH): 

^Po  'dUi  _ 'ft&l/o  /X'Zjh? 

'754-  w  r  fZ.  Oru  t- 


(2.52) 


(2.53) 


(2.54) 


(2.55) 


(2.56) 


?V; 


f?-  2%' 


(2.57) 
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Note  that  the  specific  enthalpy,  hQ,  and  the  density  may  vary 
with  both  y!  and  r,  despite  =  <$>,(_<:) . 

ORDER  S1  (FIRST)  : 

1 

It  =  -ri&W.v o)/r 

1 

(2.58) 

■ — ■ 

To  lowest  order  in  the  perturbation  quantity,  £. ,  one  may 
collect  the  following  differential  system,  written  in  convection 
form. 

1 

2S>  _u  '$£1  -  S- 

(2.59) 

[ 

f„  .  c  'gvi  .  ✓-  -jVg  _  aV. 

-5t  "  7? 

(2.60) 

i 

ft,"31**  +  P  '3lt'  1>U<  4  '&<  ^L^Z[/a 

Zb  <•<><  %  =  l  %*,»<-£  2FSy, 

(2.61) 

— 

To  3tr^  -ar  -3i5(=  R.  ’Sy.+  at 

(2.62) 

1 

Note  that  Gi  *  ^oul  here*  With  PQ(t)  imposed  externally,  (as 
expected),  this-*  system  forms  a  closure  for  $Q,  hQ,  VQ  and 

The  only  parameter  appearing  explicitly  is  the  Prandtl  number, 
which  is  of  order  unity. 

^  ~  "* 

[ 

The  associated  boundary  data  are  as  follows.  At  y^=0 
head-end  plane; 

,  the 

vQ  (0,r,t)  -  ux  (0 ,r , t)  =  0 

(2.63) 

• 

hQ  (0 ,r , t)  =  Hw  (r , t) 

(2.64) 

i 

where  Hw  is  an  arbitrary  function.  Note  that 

-  - 

fo(0,  r.t )  =  ktt)/ HJ.r.t) 

(2.65) 

t  <jC_ 

which  follows  from  the  equation  of  state,  with  p0(t)  impressed 
upon  the  layer  by  the  outer  -  flow  field.  Now  in  the  radial 
direction,  the  boundary  data  are  as  follows.  At  the  injected 

i 

-  - 

porous  wall,  r=l: 


v0{yi,l,t)  -  -Vw  (2.66) 

UiiVlflft)  =  0  (2.67) 

hQ(yi,l,t)  *  Hw(t)  (2.68) 

On  the  centerline,  r=0: 

vo(yi,0,t)  =  0  (2.69) 

"ciu^/Br  (yi,0,t)  *  0  (2.70) 

^h0/3r  (yi,0,t)  =  0  (2.71) 


Evidently,  the  second-order  derivatives  in  Eqs.  (2.60)  through 
(2.62)  are  in  the  axial  direction  (y^) ,  so  that  the  velocities 
and  thermal  enthalphy  must  match  their  outer-field  counterparts 
at  an  intermediate  region  of  common  validity. 

In  summary,  the  two  important  results  of  lower-order 
analysis  herein  are  as  follows. 

(1)  Pressure  is  uniform  within  the  layer  up  to  0{  £,1 2)  -  at 
least,  i.e.. 


Po  "  Po<fc>  '  Pi  =  Pl(fc) 

which  can  be  expected,  in  view  of  the  thin  layer  assumption 
as  well  as  the  low  velocities. 

(2)  The  axial  velocity  is  small  and  of  order  £,  ,  namely, 

u0(V;L,r,t)  =  0 

while  u^  f  0  in  general. 

The  zeroth-order  differential  system,  Eqs.  (2.59)-(2.62) , 
with  the  boundary  data,  Eqs.  (2.63)-(2.71) ,  form  a  closure. 
Obviously,  the  pressure  p_(t)  is  impressed  upon  the  layer 
externally.  Also,  external  boundary  data  is  required  for  solution 
for  the  two  momentum  equations  and  the  enthalpy  equation,  as 
expected;  this  would  enter  through  inner-outer  matching.  The 
actual  method  of  solution  will  not  be  discussed  herein. 


2.3  THE  INJECTED  SIDEWALL  LAYER 


2.3.0  THE  SINGULAR  PERTURBATION  SYSTEM 

A  porous,  injected  cylindrical  pipe  is  considered,  as  shown 
in  Fig.  3.  The  flow  region  of  interset  is  close  to  the  surface, 
where  viscous  forces  are  expected  to  be  appreciable  within  a  thin 
layer . 


For  the  neighborhood  of  r=l,  the  following  transform  is 
proposed  for  the  radial  coordinate: 


y  =  (l-r)/£ 


(2.72) 


which  magnifies  the  wall  layer,  with 

0  <  £  5  1/  <<  ! 
as  defined  earlier,  in  Eq.  (2.24).  Thus, 


(2.73) 


V3r=  -g-Way 


(2.74a) 

(2.74b) 


and 


r  =  1  -  5  y 


(2.75) 


Again  the  assumption  for  small  injection  Mach  number  is 
constrained  by: 


ocO 


(2.76) 


in  agreement  with  the  available  experimental  data. 


The  independent  variables  are  (x,y,t),  while  the  associated 
dependent  variables,  in  the  wall-layer,  are  perturbed. 


f-P.+ef.  >  v=v.-t-ev, 


V-  —VU>  +  £'Ui 


K  Shi  (2.77) 

the  following  abbreviations  are  introduced, 

£=f.Vo_,  6„=S ’,1a,  ?  }  Fj  Vo  (2.78, 


-1 


These  are  not  to  be  confused  with  their  head-end  counterparts; 
although  the  notation  is  the  same,  the  functional  dependences  are 
quite  different,  obviously. 

For  reasons  explained  fully  in  the  head-end  wall  layer 
analysis,  no  further  timescales  or  (axial)  length  scales  need  to 
be  introduced. 

2.3.1  DEPRIVATION 

For  the  perturbation  variables  of  Eqs.  (2.7 2)  —  ( 2.78) ,  the 
continuity  equation  is: 

St-fe-'-gf,)-  ^  +£*'f(Vi)  +  lfc(6o+£<&t) 

=  -((-<-6  y)(G,+eF,) 


(2.79) 


The  following  hierarchy  is  collected: 
ORDER  1/g 

—Tfefoj  =o 

E  =  fv  (x/tr) 

ORDER  6°  (ZEROTH) 

''Dg, 


—  l£‘  -_IT 
"3X  — 


TPt: 

ORDER  &1  (FIRST) 

22.  _ 

The  radial  momentum  balance  becomes: 


(2.80) 

(2.81) 

(2.82) 


(2.83) 


ire 


(&.+&,)+  dk  B(F°u'+F‘u°)]-tyc/0 

-£ %[  I? (P^fcR)J  =  -(hey)D„+  Z>f  3 


De=  £  K,  +  6  (Fa  V,r^v0)^ez(F,  «,)  V, 


> 
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n,  =  t&z(  i+e-y) 


The  hierarchy  obtained  is: 

ORDER  1/  £>3: 

—  Kut'Bfi 'of%4  =0 


Thus, 


fi>  =  fo  (*/t  ) 


ORDER  1/  &2: 

—  K^'SR/ay  =0-*'  ^  = 


(2.86) 


ORDER  1/  &  i 

—  ~c)  FbVo  /'dlsj  =  O 

ORDER  S°  (ZEROTH) : 


■Zt* 


'Sf icu 


±Vv. 

i  -21)*  . 

(2.88) 


ORDER  &  1  (FIRST) : 


75+  £(nuli-f;u.)-%[vl(Flrr,  «,)]  =  -r=oW 


(2.89) 
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The  axial  momentum  balance  is: 


££><)-*-  ^  [&0U0  i-  ^(P0t€P,) 

~  fe  ay  =“ 

+  f  0+ev)(-kX%<-£‘%>)+'^+e*$+fet%P 

~h  G^tk )%&y  (vo-+  ev.) . 

C^0V0  “♦* £(<=>oV ,  -t- <^i Vo^4- ^ (£>< V/, V<?^  ,  (2.90) 

The  associated  hierarchy  is: 

ORDER  1/  &3: 

’  K m7t&t7fK=0 


hence. 


Pc  -ft  CO 


(2.91) 

(2.92) 


ORDER  1/  €  : 


Jy(fSW.J=  0  ,2.93) 

where  we  used  FQu  =  GQv0.  This  equation  is  of  great  importance, 
as  will  be  shown  later. 


(2.93) 


ORDER  £°  ( ZEROTH) : 


(2.94) 


ORDER  fc1  (FIRST): 


+  -J^(  6,14c  +  = 

=-(s,v0+6,<)v,)-<s »v©y+-  +  ~  3  fy 

_ l  5341 

3  'SXdy 


(2.95) 
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After  similar  substitution,  the  enthalpy  equation  in  the 
wall  layer  yields  the  following  hierarchy,  for 

^  30  3  t  T°  1  }l  *  (2.96) 

ORDER  1/fc 

-T =_v  ^ 

Y°  5S  <2-97> 

Note  that  on  the  left  hand  side,  according  to  Eqs.  (2.80),  (2.87) 

fey  —  o 

(since  FQ  ^  0  in  general).  Now,  according  to  Eq.  (2.92),  pQ  = 
pQ(t),  so  that  both  sides  are  identically  zero,  as 

<4  =  4>o(t)  . 

Therefore,  Eq.  (2.97)  does  not  yield  any  new  information. 

ORDER  6°  (ZEROTH) 


If’-*-  &Wo  -  w.)  =  -  Wo  +^'5s§t 

+-Ct-i)[Uo'^  -(v,  ^  ^  )]  . 

(2.98) 

ORDER  £  1  (FIRST) 

(2.99) 


This  concludes  the  derivation  of  the  perturbed  equations  of  motion. 


2.3.2  ANALYSIS- SIDEWALL  LAYER 

The  results  of  lower-order  analysis  can  be  summarized  as 
follows.  From  Eq.  (2.81), 


f.v.  = 

while  from  Eq.  (2.87),  using  the  last  equation, 

vQ  =  vQ  (x,t)  (2.100) 

Thus  vQ  and  F  are  all  independent  of  y.  Further,  from  Eqs. 
(2.85)  and  (2.91)  clearly 

while  from  Eq.  (2.86), 

Pi  =  Pi(x,t);  (2.101) 

so  that  both  pQ  and  are  independent  of  y,  but  pG(t)  is  uniform 
within  the  entire  chamber,  as  exoected.  As  a  conseauence  of  Eq. 
(2.101), 

gk,  -t- (2.102) 

One  may  now  proceed  to  solve  Eq.  (2.93)  directly  for  uQ: 


(2.103) 

with  the  boundary  condition,  (no-slip) : 

uQ  (x,0 , t)  =  0  (2.104) 

Of  course,  the  (x,t)  dependence  of  uQ  still  remains  to  be  found. 
However,  its  dependence  upon  the  layer  coordinate,  y,  is  found 
to  be  linear;  this  result  is  certainly  not  obvious,  and  has 
several  important  implications. 


The  shear  stress  within  the  layer, 


(2.105) 
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I 

■ 


i 
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is  obviously  nonzero  in  general,  while  the  associated  azymuthal 
vorticity  within  the  layer  is  independent  of  distance  from  the 
wall  at  any  given  station  (x,t).  Even  more  striking  is  the 
vanishing  of  the  viscous  dissipation  term  at  zeroth  order: 


(2.106) 


which  leaves  in  the  zeroth-order  axial  momentum  equation  a 
balance  of  inertial  terms,  strictly,  ct.  Eq.  (2.94).  This 
explains  Dhysically  the  success  (up  to  first  order)  of  modeling 
this  family  of  injected  flows  by  assuming  rotational,  inviscid 
motions;  such  modeling  indeed  obtains  solutions  for  the  axial 
velocity  profile,  which  satisfy  the  no-slip  condition  at  the  wall 
(r=l) . 

It  further  appears  that  the  shear  stress,  Eq.  (2.105),  is 
proportional  to  the  first-order  axial  pressure  gradient,  while 
being  inversely  proportional  to  the  injected  mass  flux,  as  would 
be  expected.  Of  course,  depends  on  FQ,  and  one  expects 

their  ratio  to  be  finite  at  the  limit  as  zero  injection  is 
approached. 

In  addition  to  the  foregoing  result  for  axial  velocitv, 
the  zeroth  order  formulation  can  be  utilized  to  solve  for  the 
y-dependence  of  the  other  dependent  variables.  The  continuity 
equation  can  be  written  now  as 

^  +-  F°  =0 

From  Eqs .  (2.81)  and  (2.100)  we  know  that  F0  and  are 
independent  of  y.  Hence,  one  may  split  the  foregoing  eauation, 

'sP./3t+  FoOstr)  =  (2.107) 


=_cos<o  (2.108) 

where  CQ(x,t)  is  a  common  separation  parmeter,  with  a  range  of 
values  fully  determined  by  the  boundary  data. 

The  second  equation  yields 

|^y2/£  (2.109) 

(2.110) 
(2.111) 


where: 


FL(y,x,t)  *  B0(x,t)  +  CQ(x,t).y+ 
B0(x,t)  =  F1(0,x,t) 

,  a 

=  v: 


Similarly,  the  zeroth  order  radial  momentum  equation  is  split: 


+■  £=>v<>  =  c,  (xyt) 

lk(  ^  (y.F,  -wt  fJ)~  —  c  <,  cx  jt ) 

The  last  equation  can  be  integrated, 


(2.112) 


(2.113) 


vi(U*>0  =  B,(x;tO+-  c‘-  •y  +•  7?  ¥ 


8,(XjtO  =  v,  (Oj  x jt) 


(2.114) 


The  foregoing  results  for  and  v-^  yield  for  the  first- 
order  density: 

(^-foV.VVe  =(S0-€6lVV«+' 

-f(2Co  -C/Vc)M  a-  r^-2^  laf2- 

1  +  ^  ax  J  ^(2.115) 

Note  that:  ^(O^t)  =  (BQ  -  ^qB^Vq 

so  that  is  uniquely  defined  and  an  additional  integration 

constant  is  not  necessary.  Also,  according  to  Eqs.  (2.107)  and 

(2.112) , 


6^v./a t  =  c,_o,v„ 

Now,  since  4>£?,h.  ,  and  we  already  found  that 

■3*.  fa  =  =  3p„  (3y  =  0 


(2.116) 


then: 


hl(YrXrt)  * - gr^(y,x,t)  +  B2(x'fc) 


(2.117) 


Thus,  hi<y,x,t)  is  second  order  in  y,  like  the  appropriate 

boundary  condition  is, 


(2.118) 


B2(x,t)  =  h1{Q,x,t)  -  h0(B0  -  f0Pl)/P0 

For  the  axial  momentum,  after  dividing  through  by  F«u«,  Eg. 
(2.94)  leads  to:  °  °  ' 


JL  D&dJLo  +^02* 
•<~V4  7k  T  fb  dx 


i- 


7  'TUo 
Vc  70d 


1_  l&i  _  i_  IB  _  5^  ^AU.o 

Mo  3UJ  ^V)  fr  -3V) 


(2.119) 


Now,  according  to  the  foregoing  results: 

^'v|  =  Co/Fo 

~klk 

-  =  —  ^  /fLy 

ro  * 

Substitution  into  Eq.  (2.119)  yields,  after  some 
manipulation: 


where 

l£(*»t)s  '“•fa  ■ 

thus, 

*  (iH  -  -  BXty/f; 

and 


(2.120a) 

(2.120b) 

(2.120c) 


(2.121) 


(2.122) 


(2.123) 


uL(0,x,t)  =  0 

satisfying  the  no-slip  condition  at  the  wall.  Thus,  the 
perturbed  axial  velocity  is  third  order  in  its  y-deoendence,  and 


-27- 


the  corresponding  viscous  dissipation  term  (unlike  its  zeroth- 
order  counterpart),  does  not  vanish. 

The  zeroth  order  energy  equation  can  now  be  written  as 

^  Vo  (2.124) 

'~Mt  hr. .  —  Cl-V.Co 

<7V(raij- 

Hence,  the  right-hand  side  of  Eq.  (2.124)  becomes,  after 
substitution. 


2U» 

2* 


(2.125) 


Now  the  left  hand  side  of  Eq.  (124)  depends  only  on  t;  it 
therefore  remains  that  the  term  in  square  brackets  formally 
vanish.  Thus, 


^5U 

7SK 


TJe 

2PC 


XJcNo  —  Q  Ct ) ;  v*  f  o 


(2.126) 


One  may  turn  now  to  the  first-order  energy  equation,  which  seems 
to  yield  some  simple  and  highly  useful  rseults  even  without  full 
solution.  Equation  (2.99),  written  in  terms  of  pressure,  reads: 


+r£(p,Uc- t-fcuO-'vJijRv,  =  ->[tjv„u)+(f.V(,tp0v()] 

-  jq.  i 


JL 


(2.127) 


After  some  manipulation  one  obtains: 


H-  Uo^/TfX.  —  del'll  tfyL  —  0  (2.128) 

The  first  bracketed  term,  after  using  Eas.  (2.125)  and  (2.126), 
is  simply 

vo  ”  ^l”vo^o^  /Fo 

The  enthalpy  term,  according  to  Eas.  (2.115)  and  (2.117)  is: 


"S!h.  _  &s  _  ?l£ .  \  =w  3-rVi/^ 

~  v/*  Fo  ^  Vkk0'^)  . 


The  axial  velocity  gradient  is,  from  Eq.  (2.123): 

3tv»<  =  M  v:  H  vVi  + 

<?£>& 


(2.129) 


(2.130) 


Substitution  of  the  last  results  along  with  the  appropriate 
expression  for  Vi,  into  Eq.  (2.128)  and  collection  of  equal 
Dowers  of  y  yield's: 


JL-  ^Pl 
irpo 


+- 


-f 


Oj  Vk  ~&Jo  ,  I - 

r°l  Vo  ^  +  7»c[  vc 


Ir  M-  CoU 


~h 


4S  t  f  Kt)}^3  =0 


(2.131) 


Compatibility  with  the  foregoing  derivation  (in  which  y  and 
(x,t)  variable  separation  was  implemented),  can  be  maintained, 
provided  each  of  the  bracketed  terms  in  Eq.  (2.131)  vanishes 
identically.  The  resulting  four  compatibilitv  relations  (partial 
differential)  would  determine  the  behavior  of  the  wall  sublayer 
system  up  to  the  first  order  in  &  ,  the  small  perturbation 
quantity.  However,  a  total  of  four  undetermined  coefficients  (at 
most)  should  arise  necessarily,  to  accommodate  coupling  with  the 
outer,  inviscid  (core)  flowfield. 
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Of  particular  interest  in  the  present  analysis  is  the 
pressure, 

p(y,x,t)  *  Pd(t)+  £Pi(Xft) 

which  is  a  directly  measurable  quantity.  From  the  axial  momentum 
balance  in  perturbed  form,  cf.  Eq.  (2.90),  it  is  evident  that  the 
rotational  ("inviscid")  coreflow  can  not  sustain  a  first  order 
term  like  "bp^/dx  herein?  the  lowest-order  axial  pressure  gradient 
effect  evolves  only  at  second  order,  or  £2P2  level.  This  is 
clearly  borne  out  in  the  analyses  of  Culick,  and  others,  in  which 
the  axial  pressure  drop  is  proportional  to  M02  (Mach  number  of 
injection,  squared),  or  to  according  to  the  convention 

employed  here,  cf.  Eq.  (2.76). 

This,  however,  is  not  what  is  observed  in  the  recent  injected 
cold  flow  experiments  of  Brown,  et  al  at  CSD/UTC?  the  measured 
axial  pressure  profiles  clearly  indicate  variation  of  order  MQ  ~  £  , 
or  first  order. 


It  therefore  seems  that  the  viscous  wall  layer,  with  its 
inherent  first-order  dissipative  processes,  impresses  this  axial 
pressure  variation,  at  first  order,  over  the  entire  cross  section 
of  the  injected  channel. 

To  resolve  the  axial  variation  of  by  the  wall  layer 
formulation,  the  second  compatibility  condition  in  Eq.  (2.131) 
can  be  used,  corresponding  to  the  y-term: 


Bo  \ 

) 


-HV0+- 


Ci 


Ik  IB.  _  n 

irpo  '  O 


(2.132) 


For  the  special  case  of  uniform  (zeroth  order)  injection  at 
steady  state,  the  presence  of  a  nonzero  first-order  pressure 
perturbation  would  imply  physically  a  corresponding  nonzero 
perturbation  upon  the  mass  flux  injected,  i.e., 

Bq  =  F1(Orxvt)  /o 

as  given  by  Eq.  (2.110).  Now  at  steady  state,  although  BQ  is 
expected  _to  vary  with  p1#  we  have  assumed  for  simplicity  that 
Bq(x)  =  — Bq  =  const. 

With  the  foregoing  steady  state  assumptions.  Eg.  (2.132)  is, 
for  all  practical  purposes,  an  ordinary  differential  equation, 
for  0<x<L: 


(2.133) 
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where  "p^x)  is  the  steady  state  first-order  pressure 
perturbation;  the  coefficients  are: 


(2.134) 


C/fi o 


Note  that  at  steady  state,  according  to  Eqs.  (2.107)  and  (2.112), 
respectively,  C0=F0  and  C^=F0v0;  thus,  in  Eq.  (2.132),  C^-CQvo=0. 

The  boundary  data  are. 


dp1/dx(0)=0,  and  p1(x=L)=pL 

The  solution  is  straightforward, 

df{/JL)c  =  —  ('Ibolo',  *0 


(2.135) 

(2.136) 


(2.137) 

This  concludes  the  derivation  of  the  injected,  viscous  wall 
layer,  up  to  second  order.  Full  solutions,  namely,  matching 
between  inner  and  outer  expansions  will  not  be  attempted  herein. 
Important  insights  are  obtained  already  from  resolving  the  near¬ 
field  behavior  up  to  second  order,  in  terms  of  the  y-polynominals. 


2.3.3  DISCUSSION  OF  RESULTS 


To  facilitate  comparison  with  the  experimental  data 
reported  by  Brown,  et  al,  one  may  form  the  normalized  axial 
pressure  differential, 

Ap,  a  AR/6Pc  =  41  fe=o)-  P(*)]  = 

=  -C€+£t>  x| = 

“  *1  (2.138) 


This  axial  pressure  differential  expression  is  used  to  correlate 
the  experimental  data  of  Brown,  et  al,  as  shown  in  Fig.  4. 
Clearly,  the  measured  pressure  profile  is  correlated  very  well  by 
Eq.  (2.138),  which  is  obviously  superior  to  the  expression 
attributed  to  Culick,  shown  as  well. 


It  should  be  pointed  out  that  a  single  point  of  the  data  (x; 
p^)  has  been  utilized  to  obtain  a  scale  for  the  comparison  (this 
is  necessary,  since  no  physical  input  is  available  regardinq  the 
value  of  B0,  the  perturbed  injected  mass  flux,  necessary  for 
defining  b^,  bQ) ,  along  with  p0=F0  *  vQ  ■  1,  and  Hf  ■  1.4. 

Suppose  now  that  K._=l,  and  we  select  a  value  of  Bo=60.  (This  is 
based  on  some  trial  and  error  -  but  shows  how  the  correlation  was 
obtained  without  any  regression  analysis);  then. 


bo=l/Bo=l/60, 


b^l/TB^l/l^xSO, 
=  l/Jf1B0  =  0.014 


0.012 


Two  important  observations  are  therfore  demonstrated:  (1) 
axial  pressure  variation  to  lowest  order  is  0(£),  and  is 
governed  by  the  dissipative  wall  layer  processes,  as  shown  in  the 
rigorous  analysis  herein.  The  behavior  obtained  in  x  differs 
from  the  parabolic  pressure  droD  formula  of  Culick  [1],  and  (2) 
one  need  not  invoke  local  turbulence  generation  or  turbulence 
encroachment  upon  the  surface  to  explain  the  departure  of 
measured  Ap^  from  a  laminar  behavior. 

Another  property  of  interest  is  the  wall  friction 
coefficient,  or  dimensionless  wall  shear  stress, 

4  -  z'z  = 

where  "u*  denotes  the  mean  axial  coreflow  velocity.  Using 
dimensionless  convention  employed  herein,  along  with  the  wall 
layer  coordinate. 


(2.139) 


as  u=2x  was  used,  for  a  cylindrical  port,  and  subscript  zero 
denotes  zeroth  order  convention. 


Now,  from  Eas.  (47)  and  (77), 


where  the  first  square  root  term  is  of  order  unity.  This 
parameter  is  plotted  against  l/2x  (which  denotes  the  ratio  of 
blowing  to  mean  axial  velocity)  in  Fig.  5.  A  nearly  linear 
relationship  is  obtained,  using  the  foregoing  coefficient  values. 
In  comparison,  the  data  obtained  by  Olson  and  Eckert  (6)  is 
considered.  Ref.  6  includes  a  plot  of  the  ratio  of  (axial 
pressure  gradient)/(mean  dynamic  axial  head)  vs  v  */u0*  *  l/2x. 
This  obtains  an  almost  linear  correlation,  as  would  be  expected 
from  a  parabolic  pressure  drop.  The  slight  curvature  however, 
particularly  apparent  at  small  values  of  l/2x  <  0.01,  can  be 
followed  only  with  the  present  formulation,  not  with  any 
parabolic  pressure  profile.  Thus,  the  first  order  pressure 
distribution,  obtained  from  the  viscous  wall  layer  analysis, 
agrees  well  with  the  measured  data  of  Brown,  et  al  [8],  while  the 
associated  wall  friction  coefficient  follows  the  same  trend  as 
that  measured  by  Olson  and  Eckert  [6]. 

2.3.4  CONCLUSIONS 

A  derivation  of  the  viscous  wall  layer  regime  has  been 
presented,  pertaining  to  the  injected  flow  in  an  axial  porous 
tube,  in  simulation  of  interior  solid  propellant  rocket  flows. 

Solutions  for  the  radial  coordinate  (or  y-dependence)  of  all 
the  dependent  variables  up  to  the  second  order  have  been 
generated,  in  oolynominal  form.  The  (x,t) -dependence  is  defined 
in  terms  of  a  relatively  simple  partial  differential  system. 

Particular  results  of  the  analysis  for  the  special  case  of 
steady  state  ares  (1)  the  first  order  pressure  perturbation  was 
solved  for  and  its  axial  distribution  is  given  explicitly;  this 
term  is  entirely  due  to  the  laminar  dissipative  wall-layer 
processes,  and  (2)  the  blown  wall  friction  coefficient  was 
likewise  defined.  Both  results  correlate  well  the  available 
experimental  data.  Finally,  (3)  the  zeroth  order  axial  velocity 
distribution  within  the  layer  is  linear  radially;  thus,  to  lowest 
order,  viscous  dissipation  is  negligible  in  the  axial  momentum 
balance.  This  indicates  why  inviscid,  rotational  solutions  (such 
as  those  of  Culick  [1]  and  others,  chosen  so  as  to  satisfy  the 
no-slip  condition  at  the  wall)  are  so  successful  in  representing 
this  family  of  flows  -  up  to  first  order. 


-33- 


REFERENCES 

1.  Culick,F.E.C.,  "Rotational,  Axisymmetric  Mean  Flow  and 
Damping  of  Acoustic  Waves  in  a  Solid  Propellant  Rocket", 

AIAA  Journal,  Vol  4,  No.  8,  August  1966,  pp.  1462-1464. 

2.  Berman,  A.S.,  "Laminar  Flow  in  Channels  with  Porous  Walls", 
Journal  of  Applied  Physics,  Vol.  24,  No.  9,  September  1953,  op. 
1232-1235. 

3.  Sir  Geoffrey  Taylor,  "Fluid  Flow  in  Regions  Bound  by  Porous 
Surfaces",  Proc.  Royal  Soc.  of  London,  Ser.  A.,  Vol.  234,  No. 

1199  (1956),  pp.  456-475. 

4.  Wageman,  W.E.,  and  Guevara,  F.A.,  "Fluid  Flow  Through  a 
Porous  Channel",  The  Physics  of  Fluids,  Vol.  3,  No.  6,  November- 
December  1960,  ppT  878-881. 

5.  Dunlap,  R.,  Willoughby,  P.G.,  and  Hermsen,  R.W.,  "Flowfield 
in  the  Combustion  Chamber  of  a  Solid  Propellant  Rocket  Motor", 
AIAA  Journal,  Vol.  12,  No.  10,  October  1*974,  pp.  1440-1442. 

6.  Olson,  R.M.,  and  Eckert,  E.R.G.,  "Experimental  Studies  of 
Turbulent  Flow  in  a  Porous  Circular  Tube  with  "niform  Fluid 
Injection  Through  the  Tube  Walls",  Jour,  of  Appl.  Mech./ 
Transactions  of  the  ASME,  March  1966,  pp.  7-17. 

7.  Huesman,  K.,  and  Eckert,  E.R.G.,  "Studies  of  the  Laminar 
Flow  and  Transition  to  Turbulence  in  Porous  Tubes,  with  Uniform 
Injection  Through  the  Tube  Wall",  (Translation),  Warme  und 
Stof f uber trannunq ,  Bd.  1  (1968),  pp.  2-9. 

8.  Brown,  R.S.,  Waugh,  R.C.,  Willoughby,  P.G.,  and  Dunlap,  R., 
"Coupling  Between  Velocity  Oscillations  and  Solid  Prooellant 
Combustion",  19th  JANNAF  Combustion  Meeting,  October  1982,  CPIA 
Publication  366,  Vol.  I,  pp.  191-208. 

9.  Yagodkin,  V.I.,  "Use  of  Channels  with  Porous  Walls  for 

Studying  Flows  Which  Occur  During  Combustion  of  Solid 
Propellants",  Proc.  18th  Aeronautics  Congress,  Vol.  3,  1967,  pp. 

69-79.  (Trans.) 

10.  Yamada,  K.  Goto,  M.,  and  Ishikawa,  N.,  "Simulative  Study  on 
the  Erosive  Burning  of  Solid  Rocket  Motors",  AIAA  Journal,  Vol. 
14,  No.  9,  September  1976,  pp.  1170-1177. 

11.  Varapaev,  V.N.  and  Yagodkin,  V,I.,  "Flow  Instability  in  a 
Channel  with  Porous  Walls",  AN,  SSSR,  Mekhanika  Zhidkosti  i  Gaza, 
Vol.  4,  No.  5,  pp.  91-95,  1969.  (Trans.) 

12.  Sviridenkov,  A. A.,  and  Yagodkin,  V.I.,  "Flow  in  the  Initial 
Sections  of  Channels  with  Permeable  Walls",  lzv.  AN,  SSSR, 
Mekhanika  Zhidkosti  i  Gaza,  No.  5,  pp.  43-48,  September-October 
1976  (Moscow).  (Trans.) 


13.  Goldshtik,  M.A.,  and  Sapozhnikov,  V.A.,  "Laminar  Flow 
Stability  in  a  Mass  Force  Field",  lzv.  AN.  SSSR,  Mekhanika 
Zhidkosti  i  Gaza,  Vol.  3,  No.  5,  pp.  42-46,  1968. 

14.  Alekseev,  Yu.  N.,  and  Korotkin,  A.I.,  "Effect  of  Transverse 
Stream  Velocity  in  an  Incompressible  Boundarv  Layer  on  the 
Stability  of  the  Laminar  Flow  Regime",  Mekhanika  Zhidkosti  i 
Gaza,  Vol.  1,  No.  1,  pp.  32-36,  1966. 

15.  Flandro,  G.A.,  "Nonlinear  Time  Dependent  Combustion  of  a 
Solid  Rocket  Propellant",  Proc.  19th  JANNAF  Combustion  Meeting, 
October  1982,  CPIA  Publication  366,  Vol.  II,  pp.  111-122. 

16.  T'ien,  J.,  "Oscillatory  Burning  of  Solid  Propellants 
including  Gas  Phase  Time  Lag",  Combustion  Science  and  Technology 
Vol.  5,  1972,  pp.  47-54. 


BJ  READ-END 
BOUNDARY 
REGION 
(VISCOUS) 


9> 


INJECTED  SIDEWALL 
LAYER  (VISCOUS) 


k  «  # /  AVi*  *  *  *  j/V  °  t*^**,*  ,  O  *  •  °  «®  O  t  o  9  °o<>  Ofl  °  °  ’o  #«Bl 


COLD  FLOW 
INJECTION 


COREFLOW  CONTROL 
VOLUME  (WAVE  MOTIONS, 
MEAN  FLOW) 


CHOKED  NOZZLE 
At  =  F(t) 


FIGURE  2.1  SIMULATED  (AXI SYMMETRI C,  NONSTEADY)  ROCKET 
CHAMBER  FLOW,  SHOWING  THE  3  SPECIFIC  REGIONS 
OF  INTEREST  IN  THE  PRESENT  ANALYSIS. 
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FIGURE  2.2  SCHEMATIC  OF  THE  HEAD-END  LAYER, 

WITH  THE  STRETCHED  AXIAL  COORDINATE. 
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FIGURE  2.3 


SCHEMATIC  OF  THE  INJECTED  WALL-LAYER  REGION 
WITH  THE  PERTURBED  (EXPANDED)  COORDINATE. 
(THE  "COREFLOW"  REGIME  IS  ABOVE,  TOWARD 
THE  CENTERLINE.) 


FIGURE  2.4  DIMENSIONLESS  FIRST  ORDER  AXIAL  PRESSURE 

DISTRIBUTION,  FROM  THE  NEAR-FIELD  ANALYSIS 
HEREIN.  THE  CURRENT  RESULT  IS  PLOTTED  OVER 
THE  ORIGINAL  FIGURE  OF  CSD/UTC  (1982),  WHICH 
ALSO  SHOWS  THE  PARABOLIC  FORMULA  OF  CULICK. 


(PRESENT  STUDY) 
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FIGURE  2.5  WALL  FRICTION  COEFFICIENT ,  FROM  FIRST  ORDER 
NEAR-FIELD  ANALYSIS  HEREIN.  THE  SLIGHT 
CURVATURE  OCCURS  ALSO  IN  THE  EXPERIMENTAL 
RESULTS  OF  OLSON  AND  ECKERT  (1966);  ADJUSTMENT 
WAS  NOT  ATTEMPTED. 
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NOMENCLATURE 

=  nozzle  throat  area  and  port  exit  area,  respectively 
=  adiabatic  velocity  of  sound 
=  wall  friction  coefficient,  Eg.  (2.139) 

=  isochoric  and  isoharic  specific  heats  (J/kg-K) 

=  radial  mass  flux  (dimensionless) 

=  axial  mass  flux  (dimensionless) 

*  thermal  enthalpy,  dimensionless 

s  ratio  of  inverse  Reynolds  number  and  Mach  number 
souared,  Ecr.  (2.30),  (2.76) 

=  chamber  length 

=  Mach  number 

=  pressure 

=  Prandtl  number,  Eq.  (2.8) 

=  channel  radius 

*  injected  Reynolds  number,  Eg.  (2.7) 

=  radial  coordinate 

s  "source"-terms  in  the  eauations  of  motion  for 
coreflow,  Eqs.  (2.13)- (2.15) 

=  Strouhal  number,  injected,  (Eq.  2.34) 

=  time  (dimensionless) 

a  parameter  defining  (x,t)  -  variation  of  wall  layer 
axial  velocity  component 

=  axial  velocity,  and  mean  axial  coreflow  velocity 
respectively 

s  radial  velocity  component 

=  axial  distance  (dimensionless) 

*  radial,  magnified  wall  layer  coordinate, 
perpendicular  to  surface,  Eq.  (2.72) 

s  axial  magnified  coordinate  (head-end  boundary 
layer) ,  Eq.  (2.25) 


h 


-41- 


Greek  Symbols: 

HS  *  Cp/Cv  specific  heat  ratio 

=  difference,  increment,  Eq.  (2.138) 

=  length  scales,  Eos.  (2 . 31) - ( 2 . 34) 

=  small  perturbation  quantity,  Eqs.  (2.24);  (2.73) 
=  thermal  conductivity  of  qas  (air) ,  J/K-m-s 
=  viscosity  coefficient,  kq/m-s 
=  density 


A 
S 

£ 

A 

H- 
? 

Subscripts,  Superscripts: 


<  >o 

<  >1 
(  )* 


=  denotes  zeroth  order  (perturbation) 
=  denotes  first  order  perturbation 
=  denotes  dimensional  auantitv 


3 .  NUMERICAL  SIMULATION 


A  comprehensive  numerical  algorithm  has  been  derived  an 
implemented  for  simulation  of  the  axisymmetr ic,  nonsteady 
internal  flow  field. 

The  finite  difference  method  used  is  a  modified  MacCormack 
explicit  algorithm  (11 »  utilizing  the  original  predictor- 
corrector  scheme.  Unlike  the  original  MacCormack  algorithm, 
which  utilizes  split  time  marching  [2],  the  present  scheme  is 
unsplit  (namely,  both  radial  and  axial  space  derivatives  are 
taken  into  acount  at  each  internal  predictor-step,  within  a 
single  overall  time  increment).  This  affords  better  stability, 
particularly  near  the  walls  [3]. 

The  initial  spatial  discretization  scheme  is  shown  in  Fig. 
3.1.  Preliminary  versions  employed  uniform  radial  mesh  size,  to 
save  computation  time  in  marching  toward  steady  state  from  some 
artificial  state  described  by  the  initial  data.  Since  the  time 
marching  is  explicit,  the  smallest  spatial  increment  must 
obviously  be  used  in  the  Courant-Friedrichs-Lewy  condition, 

cAt/AR^in  <  1 

where  c=a+Vmax  is  the  dimensionless  maximal,  local  characteristic 
slope. 


The  foregoing  drawback  due  to  stability  reouirements  is 
marginal  compared  with  the  great  advantage  (at  least  in  terms  of 
the  coreflow  region  simulation),  when  compared  with  implicit 
methods  which  would  require  much  more  CPU-core  space  for  setup 
and  execution. 

The  algorithm  uses  the  dimensionless  equations  of  motion  in 
conservation  form,  as  shown  in  the  preceding  Section.  It  is 
written  in  FORTRAN  IV  and  operated  on  a  DECK  mini  computer.  A 
schematic  of  the  program  morphology  is  given  by  the  block  diagram 
of  Fig.  3.2.  A  preliminary  listing  is  provided  in  Appendix  A. 
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EXIT 

PLANE 

(TO  NOZZLE) 


X 

CENTERLINE  -  CORE  AXIS  IXX 


X  =  cAX 

typical  VALUES:  JXX  *  10,  IXX  =  25,  R0  =  1,  L  *  25. 


FIGURE  3.1  COMPUTATIONAL  MESH  FOR  AXI SYMMETRIC 
COREFLOW  SIMULATION 
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0  INTEGRATION  BY  MAC  CORMACK  METHOD 
C  TWO  DIMENSIONAL  SIMULATION  IN  R-X 

•I  V<  l,  J.  K>=RhC3=DENSITY  OF  GAS 

>:  c<2.  J,  K)*RADIAL  MOMENTUM»RHOG  *  RADIAL  VELOCITY 
<:  3.  J.  K> "AXIAL  MGMEN  T  UM=R  HOG  *  AXIAL  VELOCITY 

C  V<4.  J,  K)«T!i£RMAL  ENTHALPY“RHGG  *  ENTHALPY 

C  JXX "NUMBER  OF  RADIAL  DIVISIONS 
C  *  X  X “NUMBER  OF  AXIAL  DIVISIONS 

:  .'xmx»jx.<-i 


.  M/SEC 
,  KG/i1**3 
.  M 
,  M 

,  N/M**2 
.  KG/M-SEC 
.  JOULE /K-KC 
,  JOULE /K-M -St L 
,  M/SEC 


:  .AHA.  REYNOLDS  NUMBER,  PRANDTl  NUMBER AND  MACH  NUMBER: 
•:  :-AMA=RATIO  OF  SPECIFIC  HEAT=CP/CV 
0  -  €  0  «R  E  Yf  -'OLDS  NUMBER  *R  HOG  *  V2ERO  *  RSTAR  /  VISC 
►VN-PRANOTL  NUMBER -V 1  SC  CP  /  COND 
:  EMC  ssf'ACH  NUMBER *V.'ERO  /  SEND 
.  EI1-:*.2*MACH  N-JMBErt  f.QUARE«£MO  **  2 

:  :fl*cc!jrant-friedrichs-lewy  number 

C  NOTE:'  all  UNITS  IN  S  I  -  MKS 

:  .-ZERO-REFERENCE  injection  velocity  .  .  . 

:  -.-hog-reference  gas  density  . 

w  STAR  “MOTOR  DIAMETER < INNER )  . 

C  .STAR “MOTOR  LENGTH 
:  PSTAR*GAS  PRESSURE 

:  vise “COEFFICIENT  OF  VISCOSITY  OF  GAS 

C  CP-SPECIFIC  HEAT  OF  GAS (ISOBAR)  . 

CQND-COEFFIC IENT  OF  THERMAL.  CONDUCTIVITY  OF  GAS  . 

C  SSND “SPEED  OF  SOUND" <  GAMA*PSTR /RHOG )  ■*•*0.  5 . 


MAIN  PROGRAM 

DIMENSION  U(4,  5.  10),  RRO) 

CGMMGN/BLQCK1  /  GAMA,  CGI,  DR,  DX,  JXM1,  IXXM1.  MXM1,  JXX,  KXX,  MXX 
C0MM0N/BL0CK2/RE0,  PRN,  GAMl,  CG2,  CG3,  CG4.  CG5,  TODR,  TODX,  DR2,  DX2,  RR 
COMMON /'BL0CK3/U,  DT 

C0MM0N/BL0CK4/RSTAR,  RO,  XSTAR,  XO,  VSTAR. VO. PSTAR. PO.  RHOG,  SSND,  EMO 
*  .  CFL 

TIM* «0 

ITMAX-2 

CALL  SDATA 

OT-1  E-3 

CALL  POATA 

DO  5  L«1  .  ITMAX 

CALL  TIMIf4T<U.  TIME,  DT) 

CALL  PR  I NT ( U,  TIME) 

CONI  IMJE 


ST  CP 

END 


SUBROUTINE  TlMINTtU.  TIME,  DT> 


DIMENSION  U<4,  9.  10),  U8<4,  9.  10).  DF(4.  9.  10),  DQ(4.  9.  10*).  S(4,  9.  10). 
*  SB ( 4.  9,  10),  RR(5> 

COMMON /BL0CK1  /GAMA.  CGI,  DR,  DX,  JXMl,  KXM1,  MXM1.  JXX.  KXX.  MXX 
COMMON/ BLOC K2/RE0,  PRN,  GAM1.  CG2,  CG3,  CG4,  CG5,  TODR.  TODX.  DR2,  DX2,  RR 

2-STEP, 2- DIMENSIONAL  MAC  CORMACK  METHOD 
CALL  DFDRB(U,  DF) 

CALL  DODXBIU,  DG) 

CALL  SORCE(U, S) 

DO  600  M»1,MXX 
DO  6C0  J*2,  JXMl 
DO  6C0  K=2, KXM1 

UB(M,  J,  K)=U(M,  J,  K)-DT*DF(M,  J,  K)-DT*DG(M,  J,  K)+DT*S(M.  J,  K> 

CONI  IN'JE 


DO  60  3  M*1,MXX 
DO  4>03  K=*l.  WXX 
UBi  ri,  .'XX,  K )  =U <  fl.  JXX,  K) 

3  USiM.  1,  K)«U(M.  1.  K> 

DO  604  ri«l,MXX 

DO  *04  J*2,  JXMl 

ua < n,  j.  kxx)*U(M,  j,  kxx) 

4  UB(?1,  j,  i  >«u<n,  j,  t  > 

CALL  3N0RY <  US ) 

CAI..I.  DFDRF(UB.DF) 

CALL  000 XF (L'B,  DO) 

CALL  SORCEtUB, SB) 

DO  605  f1=l,MXX 
DO  605  J *2,  JXMl 
DO  605  K=2,  KXM1 

UU1,  J,  K )  -0.  5* ( U<  M.  J,  K)+UB(M,  J,  K)  )-DT/2.  *DF(M,  J,  K)-DT/2. 
■*  *D0(«,  J,  K)+DT/4  *(SfM,  J,  K)+SB(M,  J,  K)  ) 

5  CONTINUE 

CALL  3NDRY(U> 


TIMF.aTIf'iE+DT 

RETURN 


c 


SUBROUTINE  SOATA 


DIMENSION  U<4,  5,  10>.  RR<9> 

COMMON-  KL0CK1  / GAMA.  CGI,  L'R.  DX,  JXM1,  KXM1,  MXM1,  JXX.  KXX,  MXX 
COMMON'  BL0CK2/RE0,  PRN.  GAM1 .  C02,  CG3,  CG4,  CG9.  TODR.  TODX,  DR2.  DX2,  RR 
COMMON/ BL0CK3/U.  DT 

COMMON ' BLOC K4/R STAR,  RO.  XSTAR,  XO, VZERO. VO, PSTAR, PO, RHOG, SSND. EMO 

*  ,  CFL 

C  -R08LEM  DATA  AND  BOUNDARY  VALUES 
VZERO  n  0 
RH00=1  24 
VlSC-1.  £-9 
CP»1  16E3 
CONi>=l.  65E-2 
PSTAR=2  E5 

RSTAR-O.  09 
XSTAR-O.  29 
OAflA-1.  4 

SSNO*»  ( GAMA*PST AR  /RHOG )  #*0  9 
C  DIMENSION!  ESS  GROUP 
CFL*0  a 
OAMA-l  4 

REO-RHOG*VZERO*RSTAR / V I SC 
PRN»VISC*CP/CONO 
EMO  VZERO/SSND 
EM02a£M0**2 
C 

VOT 
PO  *  1 . 

RO-1 

XO«XSTAR/RSTAR 

C 

G AMI -GAMA- 1. 

GAi>?=GAM  1  /GAMA 
CAh3«GAMA-*EM02 
CGI =GAM2/GAM3 
CG2-4  /3.  /REO 
CG3  = 1  /REO 
CG4-GAMA/RE0/PRN 
CG9=G/\M1*GAM3/RE0 

C  MESH»9*9 

JXX-9 
k:<;<»9 
MXX =4 

uxm»jxx-i 

KXM1*KXX-1 

MXMl*MXX-l 

C 

DR»R0/JXM1 
0X»X0/KXM1 
TODR-2  *0R 
TOOX-2  *DX 
0R2-.DR-*DR 
DX2-DX-*DX 
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C  TIME  INCREMENT  CALCULATION 

C  t.  THE  CO'JRANT-FRIEDRICHS-LEWY  NUM9ER-CFL,  STATES  THAT  C*DT/DX  .  LE.  CFL 
C  WHERE  C -CHARACTER I STIC  VELOCITY,  AND  CFL  .  GT.  0  AND  .  LE.  1.  CFL  IS  AN  INP» 

C  2  DT-THE  DIMENSIONLESS  TIME  INCREMENT,  ACCORDING  TO  THE  CFL  CONDITION. 

C  ONE  HAS  TO  CHOOSE  THE  SMALLEST  (MIN.  OVER  K.  J>  ONE: 

C  DTXX-C FL*DX/CX 

C  DTRR=CFL*DR/CR 


ENX-XO/RO 
UMAX =2  >ENX 
CX-'JMAX+l.  /EMO 
CR-1.  *1.  /EMO 

dtx.<=:fl*dx/cx 

DTRR*CFL*DR/CR 
DT-AMIN1 ( DTXX,  DTRR ) 


U10  =  l. 

U20*l. 

V30  * 1 . 

U40  -G.4MA/  ( GAMA- 1 ) 
C 

DO  701  J-2, JXX 
701  RR  ( J ) -DR*  <  J- 1 > 

C  INITIAL  VALUE 


DO  700  K-l.KXX 
U< 1,  1,  K)«U10 
U(2,  1,  K)»0. 

U<  3,  1,  K)-’J10*2.  *DX*(K-1> 
U( 4,  t.  K)*'J40 

U(l,  JXX.  K)»U10 
U(2.  JXX,  K)»U20 
U< 3,  JXX,  K)»0. 

U< 4.  JXX.  K ) -U40 
CONTINUE 


700 

C  K-l 


DO  70S  J* 
U<1,  J.  1). 
U(2,  J,  1  )* 
U(3,  J.  1  )* 
U(  4,  J.  1 )  j 
CONTINUE 


>2,  JXX 
U10 
0 
0. 

'U40 


C  K-KXX 


DO  706  J-2.  JXM1 
U!l,  J,  KXX)-U10 
U<2>  J,  KXX  ) -U20 
U(3<  J»  HXX)-U10*2.  *DX*KXM1 
U!4,  J,  KXX)-U40 
706  CONTINUE 

C  INNER  POINTS 

DO  710  U-2.  JXM1 
DO  710  K-2, KXM1 
U( 1, J, K)-U10 
U(2.  J.  K)*U20 
U(3, J  . K)»U10*2  *DX*!K-1> 
U(  4,  J,  K)*U40 
7 tO  CONTINUE 

RETURN 

Ef4D 


SUBROUTINE  PDATA 
DIMENSION  U<4,  5.  10) 

CCMriON/BLOCKl/OAWA.  CGI,  DR.  DX,  JXM1.  KXM1,  MXM1.  JXX,  KXX.  MXX 


C0MM0N/BL0CK2/RE0.  PRN.  QAM1,  CQ2,  C03.  C04.  005.  TODR,  TODX.  DR2.  DX2.  RR 
C0:<f10N/BL0CK3/U.  DT 

C0I1MQN/BLQCK4/RSTAR.  RO.  XSTAR.  XO.  VZERO.  VO.  PSTAR,  PO.  RHOO.  SSND.  EMO 

*  .  CFL 

TIME-0. 

CALL  PRINTIU,  TIME) 

WRITE! 1. 1000)RSTAR,  RO.  XSTAR.  XO.  VZERO. VO. PSTAR, PO 
WRITE! 1.  1005) RHOO.  SSND.  REO,  PRN,  EMO.  GAMA,  CFL,  DR,  DX.  DT 
1000  FORMAT! 1H1. 8X.  'MOTOR  DI AMETER ! M> • ' .  F5.  2,  10X,  'RO! DIMENSIONLESS >* ', 

*  F5.  2/9X.  'MOTOR  LENGTH!M>« ',  F5.  2.  12X,  'XO! DIMENSIONLESS)- F5.  2/ 

*  9X,  'IN  JECTION  VELOCITY! M/SEC)-',  F5.  2.  '  VO! DIMENSIONLESS)- ',  F5.  2 

*  /9X»  GAS  PRESSURE < N/M**2 )-',  E9.  2.  3X,  'PO! DIMENSIONLESS)- '.  F5  2 > 

1005  FORMAT ! 9 X,  'GAS  0ENSITY!KG/M**3>- F5.  2/ 

*  9X,  'SPEED  OF  SOUND ! M/SEC  E9.  2/ 

*  9X.  'REO-',  E9.  2/9X,  'PRN- ' ,  F5.  2/9X,  'EMO- ' .  E9  2/9X,  'GAMA- ' ,  F5.  2/ 

*  9X»  'CFl»',F5  2/9X,  'DR- ',  F6.  3/9X,  'DX-  ' ,  F6.  3/9X,  'DT-',E9.  2) 

RETURN 

END 
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SUBROUTINE  BNDRY(U) 

DISSENSION  U<4.  9.  10) 

C  OMMCN  /  BLOC  K 1  /  0  AN  A  «  C01.  DR.  DX,  JXM1,  KXM1,  MXM1,  JXX.  KXX,  MXX 

DO  900  M-l.MXX 
DO  900  K-2,  KXM1 
U<M,  1.  K)«U(M,  2,  K> 

CONI INUE 

DO  909  M»2,  3 

DO  909  J-l.JXMl 

UU1,  J.  KXX  >=U(M,  J,  KXM1 ) 

CONTINUE 

RETURN 

END 


SUBROUTINE  DFDRB ( U, DFB ) 

DIMENSION  U(4»  9.  10).  DFB ( 4.  9.  10).  F(4.  9.  10) 

COMMON/ BLOC K1  /GAMA.  CGI,  DR.  DX.  JXM1,  KXM1.  MXM1.  JXX,  KXX,  MXX 

DO  ICO  M»l»  MXX 
DO  ICO  J-l.JXMl 
DO  ICO  K«2. KXM1 
F(l,  J,K)«U(2,  J.K) 

F<2.,  J.  K)«U(2.  J,  K)*U(2,  J.  K)/U<1,  J,  K)+CG1*UC4,  J.  K) 

F<3,  J,  K)*U(2.  J,  K)*U<3.  J.  K)/U(  1,  J,  K> 

F(  4,  J,  K )  =GAMA-=»U<  2,  J.  K )  *U<  4,  J,  K  > /U<  1 ,  J,  K  > 

CONI INUE 

DO  109  M-l.MXX 
DO  109  J»2.  JXM1 
DO  109  K»2, KXMt 

DFB(M. J, K)«(F(M, J, K)-F(M,  J-l,  K) >/DR 
CONTINUE 


RETURN 

END 
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SUSROUTINE  DQDXB(U, DOB) 

DIMENSION  U(4.  9.  10),  DQB(4,  9.  10).  0(4.  9.  10) 

COMMON/ BLOCK 1 /GAMA.  CQ1.  DR.  DX.  JXM1,  KXM1,  MXM1. JXX. KXX. MXX 

DO  2C0  fl-l.MXX 
DO  2C0  J*2.  JXM1 
DO  200  K»1 .  KXM1 
0(1.  J,  K)-U<3,  J,  K) 

0(2.  J,  K)»U<2.  J,  K)*U<3,  J.  K>/U(  1.  J.  K) 

0(3.  J,K)»,J(3.J,R)»U(3.J.K)/U(1.  J.  K  )+CGl*U<  4,  J,  K) 

0(4.  J,  K  > **OAMA»U (  3.  J,K)*U<4,  J.  K >  /U<  1,  J,  K ) 

CONI  I NOE 

DO  209  M*l»  MXX 
DO  209  J»2, JXM1 
DO  209  K=2, KXM1 

DOB  (M,  J,  K)**(G<M,  J.  K)-0<M,  J.  K-l  >  ) /DX 
CONI INUE 

RETURN 

END 


SUBROUTINE  SORCE(U.  S) 


DIMENSION  U< 4,  9.  10).  S(4.  9.  10).  V(3,  9,  10).  DVDR(3.  9.  10), 

*  DVDX ( 3<  9,  10),  DVDR2<3,  9,  10),  DVDX2(3,  9,  10),  DU4DR(9,  10), 

*  DU4DX ( 9.  10),  DVDRDX(2,  9.  8),  RR<9) 

CCMilCN/BLOCKl/GAMA.  C01,  DR,  DX,  JXM1,  KXM1,  MXM1,  JXX,  KXX,  MXX 
CGKMQN/BL0CK2/RE0,  PRN,  0AM1,  C02,  CQ3,  C04.  COD, TODR. TODX, DR2, DX2. RR 


C  DEFINE  V(f1,  J.  K) 

DO  300  M-l.MXMl 

DO  300  J-l,  JXX 

DO  300  K-l, KXX 

V<H,  J, K>»U<M+1,  J,  K)/U<  1,  J,  K) 

300  CONTINUE 

C  DEFINE  DVDR, DVDR2, DVDX.  DVDX2,  DVDRDX,  DU4DR,  DU4DX 
DO  309  M-l.MXMl 
DO  309  J*2,  JXMl 
DO  309  K-2, KXM1 

DVORCM,  J,  K)»(V(M,  J+t,  K)-V(M,  J-l,  K)  )  /TODR 

DVDR2(M,  J,  K)-(V(M,  J+l.  K)-2.  -*V(M,  J,  K>+V(M,  J-l,  K)  )/DR2 


DVDX (M,  J,  K)-(V<M,  J.  K+l  >-V<M,  J,  K-l )  )/TODX 

DV0X2(M,  J,  K)«(V(M.  J.  K  +  l  )-2.  #V(M,  J,  K )  +V < M,  J,  K- 1  >  > /DX2 


/-A  309 


CONI INUE 


DO  310  J»2. JXM1 
DO  310  K*2.  KXM1 

DVDRDXd.  J,  K)-(V(1,  J+l.  K+i)-V(l.  J+l#  K-l  >-V<  1.  J-1,K+1>  + 

•  Vd#  J-l.  K-l ) l/TQDR/TODX 

DVDRDX ( 2.  J,  K)*( V(2.  J+l.  K*1  >-V<2.  J+l.  K-l >-V<2.  J-l.  K+l  >♦ 

*  V<2,  J-l,  K-l) >/T0DR/T0DX 

0U4DR <J«K)*(U(4.  J+l.  K)-U<4.  J-l.  K )  > / TODR 

DU4DX( J,  K)"(U<4# J,  K+l )-U(4.  J.  K-l ) )/T0DX 
310  CONTINUE 


C  CALCULATE  S<I1.  J,  K> 

DO  319  J»2,  JXM1 
DO  319  K»2»  KXIil 
3(1.  J.  K>»-U<2,  J.  K)/RR(  J) 

S<2.  J,  K)»-U(2.  J.  K>*Vd.  J.  K )  /RR  ( J  >  +CG2/RR  ( J )  *  ( DVDR  <  1.  J,  K>- 

*  V(l, J, K)/RR(J) >+CG3*DVDX2d,  J,  K>+CG3*DVDRDX<2. J, K>/3.  + 

*  CG2*DV0R2d,  J.  K> 

S<3.  J,  K)*»-U<3.  J.  K)*V<  1,  J.  K >  /RR ( J )  +CG3/RR  <  J)  * f  DVDR  ( 2#  J,  K)  + 

*  DVDRd.  J.  K)/3.  >+CG3*DV0R2<2.  J.  K)+CG2*DVDX2<2, J. K>+CG3* 

*  DVDRDXd.  J,  K)/3 

S(4.  J.  K>«-GAMA*U<4,  J.  K)*V<  1.  J.  K> /RR  <  J>*GAMHM  V<  1.  J.  K>* 

*  DU40R<  J,  K  >+V< 2.  J.  K)#DU4DX(  J.  K)  ) +CG4*< DVDR <3.  J.  K> 

*  /RR(J)fDV0R2(3.  J, K)+DVDX2(3.  J.  K>  >+CG9*<4.  /3.  *(DVDX(2.  J,  K) 

■*  *42+ DVDR (.1,  J.  KM#2+(V(1.  J.  K)/RR<J)  >**2>+DVDR(2.  J.  K>**2 

*  ♦DVDX <  1 «  J.K) -**2-4.  *<Vd.  J.  K)/RR(J)*DVDR(1,  J.  K)+V(l.  J,  K>/ 

*  RR ( J )*DVDX ( 2.  J.  KM-DVDRd.  J.  K)*DVDX(2.  J,  K>  >  ) 

319  CONTINUE 


RETURN 

END 


SUBROUTINE  DFDRF(U,  OFF) 


DIMENSION  U<4,  9,  10),  DFF(4,  9.  10).FC4.  9.  10) 
CQMM0N/BL0CK1/GAMA,  C01.  DR.  DX.  JXMl.  KXM1,  MXM1,  JXX.  KXX. MXX 

DO  400  M»l.  MXX 
DO  400  J-2,  JXX 
DO  400  K-2,  KXM1 
F(  1,  J,  K)»U(2,  J,  K) 

F(2.  J,K>-‘J(2,  J.  K)*U<2.  J.K)/U(1,  J,  K)«-CQt*U<4,  J.K) 

F(3.  J.  K)-'J<2,  J.  K)*UC3.  J-  K)/U(  1.  J.  K) 

F(4,  J.  K)-GAMA*U<2,  J,  K  )-*U(  4,  J.  K)/U<  1.  J.  K> 

CONTINUE 

DO  409  M-l.MXX 
00  409  J-2, JXM1 
DO  409  K-2. KXM1 

DFF(M.  J,  K)»<F<M.  J*1.K)-F(M.  J.  K)  )/DR 
CONTINUE 


RETURN 

END 


SUBROUTINE  DODXFTU. DOF) 


DIMENSION  U<4.  9.  10).  D0F(4.  9.  10).  0(4.  9.  10) 

COMMON/ BLOCK  1  /GAMA.  C01.  DR.  DX.  JXM1.  KXM1.  MXM1.  JXX,  KXX.  MXX 

DO  900  M-l.MXX 
DO  900  J-2. JXMl 
DO  SCO  K-2,  KXX 
0(1.  J,  K)-U(3.  J.  K) 

0(2.  J.  K)«U(2.  J,  K)#U(3,  J,  K>  /U(  1,  J.  K> 

0(3,  J.  K)-U(3.  J.  K )*U( 3.  J,  K)/U(  1»  J»  K)+C0l*U(4,  J.  K) 

0(4,  J,  K)-GAMA*U<3,  J,  K)*U(4.  J,  K)/U<  1.  J,  K) 

CONTINUE 

DO  909  M-l.MXX 
DO  909  J-2,  JXMl 
DO  909  K-2,  KXM1 

DOF ( M,  J,  K)»(OfM,  J,  K+l  >-G(M,  J,  K)  )/DX 
CONTINUE 

RETURN 

ENO 


-55- 


C  M-4 


C  M»3 
903 


C  M*2 
902 

C  fi-1 
901 

904 


909 

810 

800 

819 

820 
929 

• 

830 

* 

839 

840 

* 

949 


SUBROUTINE  PRINT(U,  TIME) 


DIMENSION  U( 4*  9.  10) 

COMNON/BLOCKl/GAMA,  CGI.  OR.  OX.  JXM1.  KXM1,  MXM1,  JXX.  KXX.  MXX 

KM*(K.<X-l)/2+l 

DO  800  t1»l.  MXX 

IF (M. EQ.  1)  GO  TO  801 

IF(M.  EQ.  2)  GO  TO  902 

IF(M.  EQ.  3)  GO  TO  803 

WRITE  (1.840)  TIME 
WRITE! 1.  849) 

00  TO  804 

WRITE  (1.839)  TIME 
WRITE  (1.849) 

GO  TO  804 

WRITE  (1.830)  TIME 
WRITE  (1.849) 

GO  TO  804 

WRITE  (1,325)  TIME 
WRITE  (1.849) 

DO  810  K-l.KXX 

IF  (K  EQ  KM)  GO  TO  809 

WRITE  (1.819)  (U(M.  J,  K) ,  J*l«  9) 

GO  TO  810 

WRITE  (1.820)  (U(M.  J.  K>,  J-l,  5) 

CONTINUE 

CONTINUE 

FORMAT  (IX,  MX,  9(E14.  6)/) 

FORMAT  (IX,  'X'.  IX.  5(E14.  6)/) 

FORMAT  (1H1.1X.  'U(l.J.K)  :  GAS  DENSITY  VS  R  AND  X,  ',  6X, 
'TIME* '«  E10.  3//) 

FORMAT  (1H1.1X,  'U( 2.  J,  K)  :  RADIAL  MOMENTUM  VS  R  AND  X.  ',6X, 
'TIME*',  E10.  3//) 

FORMAT  (1H1.1X,  'U(3,  J,K>  :  AXIAL  MOMENTUM  VS  R  AND  X,  ',6X, 
'TIME- ',  E10.  3//) 

FORMAT  (1H1.1X,  'UC4,  J.K)  :  THERMAL  ENTHALPY  VS  R  AND  X,  ',6X, 
'TIME*',  ElO.  3//) 

FORMAT  ( IX,  ' . R . 

. '//) 

RETURN 

ENO 


lJL 


